#### Prepare figures describing the data

### Prepare environment

setwd("~/research/coffee-dataverse/src")
do.origspec <- T
source("load.R", encoding="iso-8859-1")

library(dplyr)
library(reshape2)
library(ggplot2)

### Prepare the data

df2 <- df %>% group_by(year) %>% summarize(total.produce=sum(total.produce, na.rm=T) / 1e6,
                                           total.harvest=sum(total.harvest, na.rm=T) / 1e6)
df2 <- subset(df2, year >= 1985)

### Figure showing correlation in concavity

df2$ddtot.produce <- c(NA, -df2$total.produce[1:(length(df2$total.produce)-2)] + 2*df2$total.produce[2:(length(df2$total.produce)-1)] - df2$total.produce[3:length(df2$total.produce)], NA)
df2$ddtot.harvest <- c(NA, -df2$total.harvest[1:(length(df2$total.harvest)-2)] + 2*df2$total.harvest[2:(length(df2$total.harvest)-1)] - df2$total.harvest[3:length(df2$total.harvest)], NA)
df2$ddtot.produce <- (df2$ddtot.produce - mean(df2$ddtot.produce, na.rm=T)) / sd(df2$ddtot.produce, na.rm=T)
df2$ddtot.harvest <- (df2$ddtot.harvest - mean(df2$ddtot.harvest, na.rm=T)) / sd(df2$ddtot.harvest, na.rm=T)

df3 <- melt(df2, id='year')
df3$outvar <- sapply(1:nrow(df3), function(ii) substring(df3$variable[ii], 7, 14)[[1]])
df3$facet <- sapply(1:nrow(df3), function(ii) substring(df3$variable[ii], 1, 5)[[1]])
df3$outvar.label <- ifelse(df3$outvar == "produce", "Quantity Produced (MT)", "Harvested Area (Ha)")
df3$facet.label <- ifelse(df3$facet == "total", "Levels (MT and Ha, x 1e6)", "Standardized Concavity")

ggplot(df3, aes(year, value, colour=outvar.label, linetype=outvar.label)) +
    facet_wrap(~ facet.label, ncol=1, scales='free_y')  +
    geom_line() + scale_x_continuous(expand=c(0, 0)) + theme_bw() +
    scale_colour_discrete(name=NULL) +
    scale_linetype_manual(name=NULL, values=c('dashed', 'solid')) +
    ylab(NULL) + xlab(NULL) +
    theme(legend.justification=c(1,0), legend.position=c(.9,.545))
ggsave("../figures/data/concavity.pdf", width=6.5, height=4.5)

cor(df2$ddtot.produce, df2$ddtot.harvest, use="complete")

## Produce histograms

ggplot(df, aes(1000 * gdd1000)) +
    geom_histogram() + theme_bw() + ylab("Municipality-year count") +
    xlab("Growing degree-days")
ggsave("../figures/data/preds-gdds.pdf", width=4, height=3)

ggplot(df, aes(1000 * edd1000)) +
    geom_histogram() + theme_bw() + ylab("Municipality-year count") +
    xlab("Extreme degree-days") + scale_y_log10()
ggsave("../figures/data/preds-edds.pdf", width=4, height=3)

ggplot(df, aes(tmin)) +
    geom_histogram() + theme_bw() + ylab("Municipality-year count") +
    xlab("Average minimum temperature (C)")
ggsave("../figures/data/preds-tmin.pdf", width=4, height=3)

ggplot(df, aes(prcp)) +
    geom_histogram() + theme_bw() + ylab("Municipality-year count") +
    xlab("Total season precipitation (m)")
ggsave("../figures/data/preds-prcp.pdf", width=4, height=3)

ggplot(df, aes(total.produce / total.harvest)) +
    geom_histogram() + theme_bw() + ylab("Municipality-year count") +
    xlab("Production / harvested area (MT/Ha)") + scale_x_log10()
ggsave("../figures/data/preds-yields.pdf", width=4, height=3)
